# All the code necessary to make those cool maps of voting by state, income, and ethnicity
setwd("d:/working/voucher/")
library(arm)
library(mapproj)

source("d:/FUN.R")
load("voucher2000.Rdata")

year <- 2000

threeway.fit <- glmer (data.v ~ z.inc.v*z.state.inc.v + z.inc.v*rvote04.v + (1 + z.inc.v | region.v) + (1 + z.inc.v | state.v) + (1 | inc.region.v) + (1 + z.inc.v| releth.v) + (1 | inc.releth.v) + (1 | inc.state.v) + (1 | releth.region.v) + (1 | releth.state.v) + (1 | inc.v) , family=quasibinomial(link="logit"))
display (threeway.fit)

fit <- rep (NA, n.state*n.inc*n.releth)
fit[!is.na(ybarw.v)] <- fitted (threeway.fit)

theta.hat.v <- fit

# Put the estimated cell means back as a three-way array
theta.hat.unshifted <- array (theta.hat.v, c(n.state,n.inc,n.releth), dimnames=namelist3)

theta.hat <- theta.hat.unshifted

# Print everything
pfround (prop.voters, 2)
pfround (ybar.weighted, 2)
pfround (theta.hat, 2)

# Create the totals, summing over all ethnic groups
ybar.weighted.2way <- array (NA, c(n.state,n.inc), dimnames=namelist2)
theta.hat.2way <- array (NA, c(n.state,n.inc), dimnames=namelist2)
n.eff.2way <- array (NA, c(n.state,n.inc), dimnames=namelist2)
for (i in 1:n.state){
  for (j in 1:n.inc){
    ybar.weighted.2way[i,j] <- wmean (ybar.weighted[i,j,], prop.voters[i,j,])
    theta.hat.2way[i,j] <- wmean (theta.hat[i,j,], prop.voters[i,j,])
    n.eff.2way[i,j] <- sum (n.eff[i,j,])
  }
}

# Make the maps for all voters and whites only
inclabel0 <- c ("under $20,000",
                "between $20,000 and $40,000",
                "between $40,000 and $75,000",
                "between $75,000 and $150,000",
                "over $150,000")
inclabel1 <- c("Income under $20,000", "$20-40,000", "$40-75,000", "$75-150,000", "Over $150,000" )


prop.voters.sum <- array (NA, c(n.state, n.releth))
for (i in 1:n.state){
  for (k in 1:n.releth){
    prop.voters.sum[i,k] <- sum (prop.voters[i,,k])
  }
}

# Make the map
pdf (paste ("vouchermaps", year, ".pdf", sep=""), height=14.5, width=14.5)
#====================
v.idx <- seq(1, 8)
h.idx <- c(0, seq(9, 13))
main.idx <- matrix(seq(14, 53), 8, 5)
pal.idx <- c(0,rep(54, 5))
foot.idx <- c(0, rep(55, 5))
total.mat <- rbind(c(h.idx), cbind(v.idx, main.idx), pal.idx, foot.idx)
par (mar=c(0,0,2,0), oma=c(0,0,3,0))
layout(total.mat, width=c(3, rep(4, 5)), height=c(1,rep(4,8),1.5))
#=====================
national.avg <- wmean(theta.hat,prop.voters)
for (k in 0:7){
  blankplot (c("All voters", "White\nCatholics", "White evangelical\nProtestants", "White non-evang.\nProtestants", "White other/\nno religion", "Blacks", "Hispanics", "Other races")[k+1], cex=2)
}
for(j in 1:n.inc){
  blankplot(inclabel1[j], cex=1.5)
}
for (j in 1:n.inc){
  for (k in 0:7){
    if (k==0){
      statemaps (colorblind (-(theta.hat.2way[,j]-national.avg), lo=-.25, hi=.25))
    }
    else {
      statemaps (ifelse (prop.voters[,j,k]>.01, colorblind (-(theta.hat[,j,k] - national.avg), lo=-.25, hi=.25), "white"))
    }
  }
}
par(mar=c(2,0,1,0))
pal( diverge_hcl(1001, h = c(60, 190), c = 200, l = c(10, 90)))
text(1, -1, "Yes", xpd=TRUE, cex=2)
text(0, -1, "No", xpd=TRUE, cex=2)
par(mar=c(0,0,0,0))
blankplot ("The state is left blank where a category represents less than 1% of the voters of a state", cex=1.8)
mtext("2000: Do you support school vouchers?", side=3, line=0.5, xpd=TRUE, cex=2, outer=TRUE)
#mtext (paste (year, ":  State-level support (orange) or opposition (green) on school vouchers, relative to the national average of ", round (100*national.avg), "% support", sep=""), side=3, outer=TRUE, line=1, cex=1.3)
dev.off()

# (7) Make graphs for the 50 states

for (k in 1:7){
  png (paste ("voucherplots", year, " ", namelist[[3]][k], ".png", sep=""), height=1000,width=1000)
  par (mar=c(2,2,1,0), tck=0, mgp=c(1.5,.5,0), oma=c(0,0,4,0))
  graph.dims <- c(7,7)
  par (mfrow=graph.dims)
  sort.state <- rev (order (rvote08))
  count <- 0
  for (i in sort.state){
    if (!(state.abb.long[i] %in% c("AK","HI","DC"))){
      count <- count + 1
      if (prop.voters.sum[i,k] > .1){
        plot (c(1,5), c(0,1), xlab="", ylab="", xaxt="n", yaxt="n", type="n",
              yaxs="i")
        if (count%%graph.dims[2]==1)
          axis (2, c(.02,.5,.95), c("0","50%","100%"), cex.axis=1.1)
        if (count > (48-graph.dims[2]))
          axis (1, c(1.2,3,4.8), c("poor","mid","rich"), cex.axis=1.2)
        abline (.5, 0, col="gray", lwd=.5)
        fit <- theta.hat[i,,k]
        pts <- ybar.weighted[i,,k]
        neff <- n.eff[i,,k]
        se <- sqrt (fit*(1-fit)/neff)
        text (3, .9, state.name.long[i], cex=1.3)
        for (j in 1:n.inc){
          lines (rep (j, 2), pts[j] + c(-1,1)*se[j], lwd=.5)
          points (j, min (.98, max (.02, pts[j])), pch=20, cex=1.5)
        }
        lines (1:n.inc, fit, lwd=2)
      }
      else {
        plot (c(1,5), c(0,1), xlab="", ylab="", xaxt="n", yaxt="n", bty="n", yaxs="i", type="n")
        text (3, .9, state.name.long[i], cex=1.3)
      }
    }
  }
  mtext (paste (year, ":  Raw and estimated support for school vouchers within each state among", namelist[[3]][k]), outer=TRUE, cex=1, side=3, line=1)
  dev.off()
}
